function [residual, g1, g2, g3] = model_dynamic(y, x, params, steady_state, it_)
%
% Status : Computes dynamic model for Dynare
%
% Inputs :
%   y         [#dynamic variables by 1] double    vector of endogenous variables in the order stored
%                                                 in M_.lead_lag_incidence; see the Manual
%   x         [nperiods by M_.exo_nbr] double     matrix of exogenous variables (in declaration order)
%                                                 for all simulation periods
%   steady_state  [M_.endo_nbr by 1] double       vector of steady state values
%   params    [M_.param_nbr by 1] double          vector of parameter values in declaration order
%   it_       scalar double                       time period for exogenous variables for which to evaluate the model
%
% Outputs:
%   residual  [M_.endo_nbr by 1] double    vector of residuals of the dynamic model equations in order of 
%                                          declaration of the equations.
%                                          Dynare may prepend auxiliary equations, see M_.aux_vars
%   g1        [M_.endo_nbr by #dynamic variables] double    Jacobian matrix of the dynamic model equations;
%                                                           rows: equations in order of declaration
%                                                           columns: variables in order stored in M_.lead_lag_incidence followed by the ones in M_.exo_names
%   g2        [M_.endo_nbr by (#dynamic variables)^2] double   Hessian matrix of the dynamic model equations;
%                                                              rows: equations in order of declaration
%                                                              columns: variables in order stored in M_.lead_lag_incidence followed by the ones in M_.exo_names
%   g3        [M_.endo_nbr by (#dynamic variables)^3] double   Third order derivative matrix of the dynamic model equations;
%                                                              rows: equations in order of declaration
%                                                              columns: variables in order stored in M_.lead_lag_incidence followed by the ones in M_.exo_names
%
%
% Warning : this file is generated automatically by Dynare
%           from model file (.mod)

%
% Model equations
%

residual = zeros(15, 1);
murkU__ = log(params(6)*((1-params(6))/y(11))^((1-params(6))/(params(6)+1/params(5))))+(1+params(5))/(1+params(6)*params(5))*params(11)*y(19)-(1-params(6))/(1+params(6)*params(5))*log(y(14));
sigrk_p__ = (1+params(5))/(1+params(6)*params(5))*params(12);
pi__ = 3.141592653589793;
ERkp1d__ = exp(murkU__+sigrk_p__^2/2)/2*(1+erf((y(18)-murkU__-sigrk_p__^2)/(sigrk_p__*1.414213562373095)))+y(10)*(1-params(7));
ERkp1u__ = 0.5*exp(murkU__+sigrk_p__^2/2)*(1+erf((murkU__+sigrk_p__^2-y(18))/(sigrk_p__*1.414213562373095)))+(1-params(7))*(1-y(10));
T68 = 1+1/params(5);
T69 = y(21)^T68;
T72 = y(20)-y(11)*T69/T68;
T75 = y(23)^T68;
T78 = y(22)-y(11)*T75/T68;
T79 = T72/T78;
T90 = 1+ERkp1d__*y(16)/(y(17)*(y(16)-1))-y(10)*params(8)-y(10);
T100 = exp(y(19))*y(1)^params(6);
T101 = y(21)^(1-params(6));
T102 = T100*T101;
T123 = params(7)+y(17)*(1-1/y(16))*(1+params(8)*(1-y(13)))-1;
T138 = params(8)*y(17)^2;
T139 = (1-y(13))*T138;
T142 = y(16)^2;
T143 = (y(16)-1)*(1+params(8)*(1-y(13)))*T139/T142;
T149 = T143/T123/(sigrk_p__*sqrt(2*pi__));
T154 = exp((-0.5)*((y(18)-murkU__)/sigrk_p__)^2);
lhs =y(10);
rhs =0.5*(1+erf((y(18)-murkU__)/(sigrk_p__*1.414213562373095)));
residual(1)= lhs-rhs;
lhs =1/(params(4)*y(17));
rhs =T79*T90;
residual(2)= lhs-rhs;
lhs =y(7);
rhs =T69*y(11)/(1-params(6));
residual(3)= lhs-rhs;
lhs =y(7);
rhs =T102;
residual(4)= lhs-rhs;
lhs =y(14);
rhs =(1-params(7))*y(1)+y(8);
residual(5)= lhs-rhs;
lhs =y(7);
rhs =y(20)+y(8);
residual(6)= lhs-rhs;
lhs =y(14);
rhs =y(16)*y(15);
residual(7)= lhs-rhs;
lhs =exp(y(18));
rhs =T123;
residual(8)= lhs-rhs;
lhs =y(17)*(1-y(10));
rhs =ERkp1u__-T149*T154;
residual(9)= lhs-rhs;
lhs =y(15);
rhs =params(10)*(y(2)+params(9)*(y(1)*(1+params(6)*y(7)/y(1)-params(7))-y(1)*y(4)*(1-1/y(3))-y(2)))+y(12)-x(it_, 2);
residual(10)= lhs-rhs;
lhs =y(19);
rhs =params(11)*y(6)+params(12)*x(it_, 1);
residual(11)= lhs-rhs;
lhs =y(9);
rhs =params(6)*y(7)/y(1)/exp(y(5));
residual(12)= lhs-rhs;
lhs =y(11);
rhs =(steady_state(5));
residual(13)= lhs-rhs;
lhs =y(12);
rhs =(steady_state(6));
residual(14)= lhs-rhs;
lhs =y(13);
rhs =(steady_state(7));
residual(15)= lhs-rhs;
if nargout >= 2,
  g1 = zeros(15, 25);

  %
  % Jacobian matrix
  %

T15 = (1-params(6))/(params(6)+1/params(5));
T24 = (1+params(5))/(1+params(6)*params(5))*params(11);
T48 = sigrk_p__^2;
T51 = exp(murkU__+T48/2);
T52 = T51/2;
T56 = 1+erf((y(18)-murkU__-T48)/(sigrk_p__*1.414213562373095));
T130 = 1+erf((murkU__+T48-y(18))/(sigrk_p__*1.414213562373095));
T221 = params(6)*(-(1-params(6)))/(y(11)*y(11))*getPowerDeriv((1-params(6))/y(11),T15,1)/(params(6)*((1-params(6))/y(11))^T15);
T228 = 2/(exp(((y(18)-murkU__)/(sigrk_p__*1.414213562373095))^2)*1.772453850905516);
T246 = 2/(1.772453850905516*exp(((y(18)-murkU__-T48)/(sigrk_p__*1.414213562373095))^2));
T265 = 2/(1.772453850905516*exp(((murkU__+T48-y(18))/(sigrk_p__*1.414213562373095))^2));
T271 = 2*(y(18)-murkU__)/sigrk_p__;
T311 = (1-params(6))/(1+params(6)*params(5))*1/y(14);
T437 = (-T24)/(sigrk_p__*1.414213562373095);
T472 = getPowerDeriv(y(21),T68,1);
  g1(1,10)=1;
  g1(1,11)=(-(0.5*(-T221)/(sigrk_p__*1.414213562373095)*T228));
  g1(1,14)=(-(0.5*T228*T311/(sigrk_p__*1.414213562373095)));
  g1(1,18)=(-(0.5*T228*1/(sigrk_p__*1.414213562373095)));
  g1(1,19)=(-(0.5*T228*T437));
  g1(2,10)=(-(T79*((1-params(7))*y(16)/(y(17)*(y(16)-1))-params(8)-1)));
  g1(2,11)=(-(T90*(T78*(-(T69/T68))-T72*(-(T75/T68)))/(T78*T78)+T79*y(16)*(T56*T51*T221/2+T52*(-T221)/(sigrk_p__*1.414213562373095)*T246)/(y(17)*(y(16)-1))));
  g1(2,14)=(-(T79*y(16)*(T56*T51*(-T311)/2+T52*T246*T311/(sigrk_p__*1.414213562373095))/(y(17)*(y(16)-1))));
  g1(2,16)=(-(T79*(ERkp1d__*y(17)*(y(16)-1)-y(17)*ERkp1d__*y(16))/(y(17)*(y(16)-1)*y(17)*(y(16)-1))));
  g1(2,17)=(-params(4))/(params(4)*y(17)*params(4)*y(17))-T79*(-(ERkp1d__*y(16)*(y(16)-1)))/(y(17)*(y(16)-1)*y(17)*(y(16)-1));
  g1(2,18)=(-(T79*y(16)*T52*T246*1/(sigrk_p__*1.414213562373095)/(y(17)*(y(16)-1))));
  g1(2,19)=(-(T79*y(16)*(T56*T24*T51/2+T52*T246*T437)/(y(17)*(y(16)-1))));
  g1(2,20)=(-(T90*1/T78));
  g1(2,22)=(-(T90*(-T72)/(T78*T78)));
  g1(2,21)=(-(T90*(-(y(11)*T472/T68))/T78));
  g1(2,23)=(-(T90*(-(T72*(-(y(11)*getPowerDeriv(y(23),T68,1)/T68))))/(T78*T78)));
  g1(3,7)=1;
  g1(3,11)=(-(T69*1/(1-params(6))));
  g1(3,21)=(-(y(11)/(1-params(6))*T472));
  g1(4,7)=1;
  g1(4,1)=(-(T101*exp(y(19))*getPowerDeriv(y(1),params(6),1)));
  g1(4,19)=(-T102);
  g1(4,21)=(-(T100*getPowerDeriv(y(21),1-params(6),1)));
  g1(5,8)=(-1);
  g1(5,1)=(-(1-params(7)));
  g1(5,14)=1;
  g1(6,7)=1;
  g1(6,8)=(-1);
  g1(6,20)=(-1);
  g1(7,14)=1;
  g1(7,15)=(-y(16));
  g1(7,16)=(-y(15));
  g1(8,13)=(-(y(17)*(1-1/y(16))*(-params(8))));
  g1(8,16)=(-((1+params(8)*(1-y(13)))*y(17)*(-((-1)/(y(16)*y(16))))));
  g1(8,17)=(-((1-1/y(16))*(1+params(8)*(1-y(13)))));
  g1(8,18)=exp(y(18));
  g1(9,10)=(-y(17))-(-(1-params(7)));
  g1(9,11)=(-(T130*0.5*T51*T221+0.5*T51*T221/(sigrk_p__*1.414213562373095)*T265-T149*T154*(-0.5)*(-T221)/sigrk_p__*T271));
  g1(9,13)=T154*(T123*(y(16)-1)*(T139*(-params(8))+(1+params(8)*(1-y(13)))*(-T138))/T142-T143*y(17)*(1-1/y(16))*(-params(8)))/(T123*T123)/(sigrk_p__*sqrt(2*pi__));
  g1(9,14)=(-(T130*0.5*T51*(-T311)+0.5*T51*T265*(-T311)/(sigrk_p__*1.414213562373095)-T149*T154*(-0.5)*T271*T311/sigrk_p__));
  g1(9,16)=T154*(T123*((1+params(8)*(1-y(13)))*T139*T142-(y(16)-1)*(1+params(8)*(1-y(13)))*T139*2*y(16))/(T142*T142)-T143*(1+params(8)*(1-y(13)))*y(17)*(-((-1)/(y(16)*y(16)))))/(T123*T123)/(sigrk_p__*sqrt(2*pi__));
  g1(9,17)=1-y(10)-(-(T154*(T123*(y(16)-1)*(1+params(8)*(1-y(13)))*(1-y(13))*params(8)*2*y(17)/T142-T143*(1-1/y(16))*(1+params(8)*(1-y(13))))/(T123*T123)/(sigrk_p__*sqrt(2*pi__))));
  g1(9,18)=(-(0.5*T51*T265*(-1)/(sigrk_p__*1.414213562373095)-T149*T154*(-0.5)*T271*1/sigrk_p__));
  g1(9,19)=(-(T130*0.5*T24*T51+0.5*T51*T265*T24/(sigrk_p__*1.414213562373095)-T149*T154*(-0.5)*T271*(-T24)/sigrk_p__));
  g1(10,7)=(-(params(10)*params(9)*y(1)*params(6)/y(1)));
  g1(10,12)=(-1);
  g1(10,1)=(-(params(10)*params(9)*(1+params(6)*y(7)/y(1)-params(7)+y(1)*(-(params(6)*y(7)))/(y(1)*y(1))-y(4)*(1-1/y(3)))));
  g1(10,2)=(-(params(10)*(1-params(9))));
  g1(10,15)=1;
  g1(10,3)=(-(params(10)*params(9)*(-(y(1)*y(4)*(-((-1)/(y(3)*y(3))))))));
  g1(10,4)=(-(params(10)*params(9)*(-(y(1)*(1-1/y(3))))));
  g1(10,25)=1;
  g1(11,6)=(-params(11));
  g1(11,19)=1;
  g1(11,24)=(-params(12));
  g1(12,7)=(-(params(6)/y(1)/exp(y(5))));
  g1(12,9)=1;
  g1(12,1)=(-((-(params(6)*y(7)))/(y(1)*y(1))/exp(y(5))));
  g1(12,5)=(-((-(params(6)*y(7)/y(1)*exp(y(5))))/(exp(y(5))*exp(y(5)))));
  g1(13,11)=1;
  g1(14,12)=1;
  g1(15,13)=1;

if nargout >= 3,
  %
  % Hessian matrix
  %

  g2 = sparse([],[],[],15,625);
if nargout >= 4,
  %
  % Third order derivatives
  %

  g3 = sparse([],[],[],15,15625);
end
end
end
end
